Dynamical criticality of spin-shear coupling in van der Waals antiferromagnets

The interplay between a multitude of electronic, spin, and lattice degrees of freedom underlies the complex phase diagrams of quantum materials. Layer stacking in van der Waals (vdW) heterostructures is responsible for exotic electronic and magnetic properties, which inspires stacking control of two-dimensional magnetism. Beyond the interplay between stacking order and interlayer magnetism, we discover a spin-shear coupling mechanism in which a subtle shear of the atomic layers can have a profound effect on the intralayer magnetic order in a family of vdW antiferromagnets. Using time-resolved X-ray diffraction and optical linear dichroism measurements, interlayer shear is identified as the primary structural degree of freedom that couples with magnetic order. The recovery times of both shear and magnetic order upon optical excitation diverge at the magnetic ordering temperature with the same critical exponent. The time-dependent Ginzburg-Landau theory shows that this concurrent critical slowing down arises from a linear coupling of the interlayer shear to the magnetic order, which is dictated by the broken mirror symmetry intrinsic to the monoclinic stacking. Our results highlight the importance of interlayer shear in ultrafast control of magnetic order via spin-mechanical coupling.


Supplementary Figure 2:
Calibrating the peak shift measured by the center of mass on the detector in reciprocal unit. a Schematic of XRD reflection geometry. b, Schematic of the Ewald sphere and Bragg peaks in reciprocal space. Green and red disk: finite size Bragg peaks before and after laser excitation. Red dashed circle and red lines: Bragg condition for t > 0. c, The Bragg peak shifts as measured by the center of mass for 002 and 2 ̅ 02 peaks in FePS3 are compared with the results obtained by the L scans. An example of L scan at 95 K is shown in Supplementary percentage change of c* as a function of temperature. The red lines represent the linear thermal expansion at temperatures below and above TN. b, Simulation of the thermal relaxation time. The gray line is the specific heat adapted from Ref. 1 . The green data points are the relaxation time measured by c* change of 002 peak. The blue curve is the simulation results based on a 1D thermal transport model as described in Supplementary Note 2. The error bars represent the standard deviation of exponential fitting. c, Thickness-dependent thermal relaxation time measured at 100 K. The red line is a linear fit of the data. The error bars represent the standard deviation of exponential fitting. Figure 5: Calculation of monoclinic angle change (Δ ) and its dependence on absorbed energy density. a, Changes of  based on 002 and 2 ̅ 02 peak shift. b,  decrease due to an interlayer expansion (see Supplementary Note 3). c, Maximal Δ at 100 K and 135 K as a function of absorbed energy density. The green dashed line indicates the threshold fluence at 100 K, as determined from the crossing of the linear trends (thin red lines) extrapolated from the lowand high-fluence limit. Figure 6: Fitting of structural dynamics of FePS3 by exponential decay functions. a & b, Single-exponential decay fitting of 002 and 02 ̅ 2 peak at 70 K, respectively. c, Biexponential decay fitting of 2 ̅ 02 peak at 70 K. The green curve corresponds to the relaxation of interlayer expansion measured by the 002 Bragg peak while red curve corresponds to the relaxation of the  angle change. Purple curve is the sum of green and red curves. d, Comparison of the fitted amplitudes for thermal expansion (green curve) and interlayer shear (red curve) measured by 002 and 2 ̅ 02 peak as a function of the temperature. Figure 7: Comparison of Δ LD and Δ dynamics in MPS3 (M = Fe, Ni, Mn). ac, Crystal and magnetic structures of FePS3, NiPS3, and MnPS3, respectively. d-f, Δ LD dynamics in FePS3, NiPS3, and MnPS3, respectively. g-i, Δ dynamics in FePS3, NiPS3, and MnPS3, respectively. The Δ LD and Δ dynamics are plotted in the same scale for the three compounds for comparison. The sample thicknesses are the following: FePS3: 60 nm (d) and 1300 nm (g). NiPS3: 40 nm (e) and 196 nm (h). MnPS3: 50 nm (f) and 1100 nm (i). The sample thicknesses were chosen for optimizing the signal-to-noise ratios of the corresponding x-ray and optical probes. It is an extrinsic factor with a linear relation to the absolute delay time. But the thickness does not affect the scaling exponents as shown in Fig.3. The pump fluence for all the OLD measurements was ~1.0 mJ•cm -2 (400 nm, 1 kHz repetition rate). The pump fluences for the XRD measurements was set such that the interlayer thermal expansion was similar (Δc*/c* ≈ 0.05%) for all three compounds.

Supplementary Note 1: Comparison of peak center of mass shift and HKL scans
The time-resolved XRD measurements were performed in such a scheme that the sample and the detector were fixed while the diffraction peak was recorded on an area detector as a function of time delays. We show below that the change of the center of mass on the detector at the fixed incident angle in our experimental setup can be used to measure the Bragg peak shift in reciprocal space.
As shown in Supplementary Figure 2a, before time zero, the Bragg condition is satisfied, i.e., = = , where α, β, are the incident angle, exit angle, and Bragg angle, respectively. After laser excitation, the change of atomic plane spacing results in a change of the Bragg condition. As shown in Supplementary Figure 2b, the green and blue disks represent the finite-size Bragg peak. Δθ is the measured angular peak shift on the detector from point A to B. ∠AOC = 2 is the 2θ angle before laser excitation. ∠BOC is the 2θ angle after laser excitation. After laser excitation, ∠BOʹC is the theoretical 2θ angle satisfying the Bragg condition, which can be calculated based on the peak shift on the detector. Based on trigonometry, one obtains In the experiment, for 002 peak is 10.15°. The measured peak shift on the detector at +5 ns corresponds to |Δ2θ| = 0.01°. The theoretical angular shift |Δ2θ| = 2 -∠BOʹC = 0.0103°, 3% larger than the measured peak shift of 0.01°. To confirm that the systematic error of measuring Bragg peak shift using diffraction peak center of mass change is negligible, the extracted shifts of 002 and 2 ̅ 02 peak were compared with the scans along the miller index L direction (L scan) using SPEC from Certified Scientific Software (Supplementary Figure 2c). As the figure shows, for both 002 and 2 ̅ 02 peaks, the peak center of mass shift is consistent with the L scan.

Supplementary Note 2: Estimation of transient temperature rise and simulation of thermal relaxation
The sample temperature rise upon laser excitation was estimated based on the interlayer lattice expansion. As shown in Supplementary Figure 4a, the measured change of c * (interlayer spacing = 2 * ) follows a linear dependence on the temperature above and below TN. The kink at TN indicates a structural phase transition. The interlayer thermal expansion coefficient was estimated to be 1.818 × 10 −5 K −1 . In the pump-probe measurements, at a static temperature of 80 K, the measured Bragg peak shift Δc * /c * = -0.055% corresponds to a transient temperature increase of 30 K, which corresponds to an incident fluence of 5.7 mJ•cm -2 on the 1.3 μm-thick sample.
A one-dimensional (1D) thermal transport model was employed to simulate the cooling of the sample including the latent heat during the phase transition, following the method detailed in Ref. 2 . Using the following parameters: sample interlayer thermal conductivity 3 = 0.85 W • m −1 • K −1 , sample-substrate interfacial thermal conductivity G=0.17 W • m −1 • K −1 , the measured specific heat curve 1 , the simulation results were summarized in Supplementary Figure 4b and captured the key features of experimental data. The thermal dissipation time gradually increases with sample temperature and peaks at 110 K, rather than TN. The linear dependence of thermal relaxation time on sample thickness was also reproduced (Supplementary Figure 4c).
The laser-induced  decrease of 0.1° is consistent with the temperature-dependent XRD results, as shown in Supplementary Figure 9.
The change due to interlayer expansion was also observed experimentally. As shown in Supplementary Figure 5c, at the sample temperature of 135 K, Δ increases linearly with the absorbed energy shown by the black curve. At the sample temperature of 100 K, the slope of  change as a function of absorbed energy (the green curve) is 10 times larger than that measured at 135 K, corresponding to the ~10 times difference between total  change and the interlayer expansion-induced  change. Above the threshold of the absorbed energy density 3.3 × 10 4 mJ • cm −3 , the slope of the green curve matches the black curve measured as 135 K, suggesting that there is no longer interlayer shear and the  change purely comes from interlayer expansion.

Supplementary Note 4: Fitting of the time-dependent phenomena
The relaxation time of XRD and OLD dynamics were extracted based on exponential function fitting. A single exponential decay function was used for fitting the time-dependent relaxation measured by 002 and 02 ̅ 2 Bragg peak, as well as OLD, while a bi-exponential decay function was used to fit the dynamics measured by 2 ̅ 02 Bragg peak.
The single exponential decay function is: TN, the rise time of 2 ̅ 02 peak is longer than 0.1 ns so the rise time τr was also a fitting parameter, which does not affect the fitting results of the relaxation time constant A representative fitting result is shown in Supplementary Figure 6c. Below TN, the interlayer shear has a much larger amplitude than thermal expansion (Supplementary Figure 6d) and dominates the relaxation process.

Supplementary Note 5: Fitting of the critical exponents
The divergence in of both and LD as a function of temperature were fitted by a powerlaw function: = |1 − | − , where A is the amplitude, z and ν are the dynamical exponent and critical exponent of correlation length, respectively, T is the sample temperature. To better determine the sample temperature for each measurement, we first fit the data in a linear scale with three fitting parameters, , TN, and zv, which yield an initial value of TN. We noted that this value has an offset of ±2 K with respect to TN, which is attributed to the systematic error in the temperature measurements under laser illumination at different experimental conditions. Then, for better fitting of zv, the value of TN was fixed in the fitting procedure using the equation in the log scale: log = − log |1 − | + log , with only two fitting parameters: and A. The fitted zν values were summarized in Supplementary Table 1 for samples with various thicknesses, with a representative fit shown in Figure 8a. The width of the temperature range in which the slowing down is observed can be determined by the temperature at which the critical fluctuation is larger than the thermal fluctuation. This region spans from 105 K to 117 K (Fig. 2e) as the recovery of lattice becomes much longer than the thermal recovery. The width of this region is wider than those observed such as in ultrathin Fe films 4 due to weak vdW interaction and large spin fluctuations in vdW magnets.  Figure 8b & c). The fitted value was ν = 0.44 ± 0.11, close to the theoretical value of ν = 0.5 (Ref. 6 ).

Measurement
As shown in Supplementary Figure 7, the shear amplitude and the magnetic order measured by OLD are larger in FePS3 than that in NiPS3. A possible explanation is the different electron configurations of the corresponding ions. For example, the Fe 2+ ion has d 6 configuration, which leaves the minority t2g manifold partially filled, but the Ni 2+ ion has d 8 configuration, so the minority t2g manifold is fully filled, which quenches the angular momentum and in general leads to weaker spin-lattice coupling compared to partially filled d-orbitals. In addition, the strength of critical slowing down is also different: the recovery time of the shear motion in FePS3 increased by two orders of magnitude and only increased by one order of magnitude in NiPS3. This is likely due to the different spin configurations in these two compounds. The zigzag spin order is Ising type and 2D XY type in FePS3 and NiPS3 respectively, which have distinct scaling exponents and universality classes.

Supplementary Note 6: Temperature-dependent single-crystal XRD for FePS3
We performed single-crystal XRD as a function of temperature at the 7ID-C beamline of the Advanced Photon Source. During the measurement, the orientation matrix was optimized and fixed at 80 K. The H, K, L values of multiple Bragg peaks (002, 004, 2 ̅ 04, 4 ̅ 24) were recorded at each temperature and used to calculate the a, b, c, and  as summarized in Supplementary Figure   9.
We note that previous reports of the structural transition in FePS3 used the powder XRD method 7,8 did not report monoclinic angle  change as a function of temperature. To cross check our single-crystal structure refinement results, we also performed high-resolution powder XRD at 11BM beamline at the Advanced Photon Source. The powder XRD results yield a distinct  change of 0.02° across TN. This is a factor of 5 smaller than the single-crystal result. We ascribe this discrepancy to the sample preparation for powder diffraction, in which the samples were grinded into powder and attached to the Kapton tubes with vacuum grease. This process likely introduced strains in the sample, resulting in a smaller  change in powder XRD data.